hide
# Set ggplot theme for figures
theme_set(theme_bw() +
theme(
# Titles
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Axes
axis.line.y = element_line(),
axis.ticks.y = element_line(),
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.text = element_text(size = 11),
# Facets
strip.text = element_text(color = "white",
face = "bold",
size = 11),
strip.background = element_rect(
color = "white",
fill = "#00736D",
size = 1,
linetype="solid"))
)
# Building a GMRI theme based on Wall street Journal and NYTimes theme
# base settings from {ggthemes}
theme_gmri <- function(base_size = 10,
bg_color = "lightblue",
base_family = "sans",
title_family = "sans",
facet_color = "teal") {
# Color from gmRi palette, sets background color
#colorhex <- gmRi::gmri_cols()[bg_color]
facet_hex <- gmri_cols()[facet_color]
# Set up theme
theme_foundation(
base_size = base_size,
base_family = base_family) +
theme(
# Major Elements
line = element_line(linetype = 1, colour = "black"),
rect = element_rect(fill = "transparent",
linetype = 0,
colour = NA),
text = element_text(colour = "black"),
title = element_text(family = title_family, size = 12),
# Axis elements
axis.text.x = element_text(colour = NULL),
axis.text.y = element_text(colour = NULL),
axis.ticks = element_line(colour = NULL),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = NULL),
axis.line = element_line(),
axis.line.y = element_blank(),
axis.text = element_text(size = 11),
# Legend Elements
legend.background = element_rect(),
legend.position = "top",
legend.direction = "horizontal",
legend.box = "vertical",
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Panel/Grid Setup
panel.grid = element_line(colour = NULL, linetype = 3, color = "gray70"),
panel.grid.major = element_line(colour = "black"),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
# Title and Caption Details
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
plot.margin = unit(c(1, 1, 1, 1), "lines"),
# Facet Details
strip.text = element_text(color = "white", face = "bold", size = 11),
strip.background = element_rect(
color = "white",
fill = facet_hex,
size = 1,
linetype="solid"))
}
# Set theme up for maps
map_theme <- list(
theme(
panel.border = element_rect(color = "black", fill = NA),
plot.background = element_rect(color = "transparent", fill = "transparent"),
line = element_blank(),
axis.title.x = element_blank(), # turn off titles
axis.title.y = element_blank(),
#legend.position = "bottom",
legend.title.align = 0.5))
2021 Sea Surface Temperature
An Exceptional Year
2021 was an exceptionally warm year for the Gulf of Maine. It had the warmest fall on record for the area, and the second warmest summer. Much of the year the region experienced what would be considered “heatwave conditions”.
hide
# Load the timeseries
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)
# Clean up the data - add labels
region_timeseries <- region_timeseries %>%
distinct(time, .keep_all = T) %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
#Set up correct year for winters, they carry across into next year
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
season_yr = ifelse( (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(
temperature_timeseries = region_timeseries,
threshold = 90,
clim_ref_period = c("1982-01-01", "2011-12-31"))
# Format dates and seasons and heatwave outcomes
region_hw <- region_hw %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
#Set up correct year for winters, they carry across into next year
season_yr = ifelse(
(season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Set up heatwave status outcomes for fancy table features:
region_hw <- region_hw %>%
mutate(hw_outcomes = case_when(
mhw_event == TRUE ~ 1,
mcs_event == TRUE ~ 0,
TRUE ~ 0.5))
Departure from Historic Conditions
Thanks to some long-running temperature monitoring programs we are able to track monthly sea surface temperatures as far back as 1850. The figure below shows how the region’s temperatures have fluctuated as measured by the Extended Reconstructed Sea Surface Temperature (ERSST) record.
Once we hit 1982 the availability of a higher resolution temperature resource becomes available in the form of the Optimum Interpolation Sea Surface Temperature (OISST) record. The ERSST data is of a coarser resolution (monthly measurements at a 0.5 x 0.5 degree grid) and does not capture the warmer inshore dynamics (a bias to show colder temperature), but is within half of a degree of the more modern measurements.
hide
# Pathto ERSST on box
ersst_path <- cs_path("okn", "ersst")
# Load temps and anoms (1982-2011 climatology)
ersst_temp <- stack(str_c(ersst_path, "ERSSTv5_merged.nc"))
ersst_anom <- stack(str_c(ersst_path, "ERSSTv5_anom.nc"))
# Crop GoM
# function to mask
mask_nc <- function(ras_obj, mask_shape){
# Get stats from ranks
m1 <- raster::mask(rotate(ras_obj), mask_shape)
m1 <- raster::crop(m1, mask_shape)
return(m1)}
# Load the shape, mask ERSST with it
gom_shape <- read_sf(poly_path)
gom_ersst <- mask_nc(ras_obj = ersst_temp, mask_shape = gom_shape)
gom_ersst_anom <- mask_nc(ras_obj = ersst_anom, mask_shape = gom_shape)
# Clean them up - temp
ersst_temp_tl <- as.data.frame(cellStats(gom_ersst, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# anomalies
ersst_anom_tl <- as.data.frame(cellStats(gom_ersst_anom, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst_anom")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# Join them together
ersst_tl <- left_join(ersst_temp_tl, ersst_anom_tl) %>%
mutate(yr = as.numeric(yr))
# Get Yearly Average
ersst_yr <- ersst_tl %>%
group_by(yr) %>%
summarise(sst = mean(sst),
sst_anom = mean(sst_anom),
.groups = "drop") %>%
mutate(yr = as.numeric(yr)) %>%
filter(yr <= 2019)
# Same for OISST
oisst_yr <- region_hw %>%
filter(between(yr, 1982, 2021)) %>%
group_by(yr) %>%
summarise(sst = mean(sst), .groups = "drop")
# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))
# Plot them
ggplot() +
geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"),
size = 1) +
geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"),
size = 1) +
scale_color_gmri() +
theme_gmri() +
labs(x = "", y = "Sea Surface Temperature",
color = "Temperature Data Record",
title = "Long-Term Temperature Record for the Gulf of Maine",
subtitle = "Current temperatures above 150-year highs") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
facet_zoom(xlim = c(1982, 2021)) +
theme(legend.position = "bottom")
Relationship to Climate Patterns
To better get a sense of what kind of year 2021 was it is helpful to relate it to what we would expect from similar years.
ENSO Index
One large feature in global weather patterns is the El Niño Southern Oscillation (ENSO). While ENSO and the southern oscillation index are derived from temperature/pressure differences in the Pacific Ocean, they influence atmospheric behavior coming across the United States which has the potential to influence SST patterns.
hide
# Source: https://www.ncdc.noaa.gov/teleconnections/enso/soi
# Does not work:
# library(rsoi)
# clim_idx <- download_soi()
# Manually copied from link above
soi_std <- tibble::tribble(#####
~"YEAR", ~"JAN", ~"FEB", ~"MAR", ~"APR", ~"MAY", ~"JUN", ~"JUL", ~"AUG", ~"SEP", ~"OCT", ~"NOV", ~"DEC",
1951, 1.5, 0.9, -0.1, -0.3, -0.7, 0.2, -1.0, -0.2, -1.1, -1.0, -0.8, -0.7,
1952, -0.9, -0.6, 0.5, -0.2, 0.8, 0.7, 0.5, 0.1, -0.2, 0.4, 0.0, -1.2,
1953, 0.3, -0.5, -0.2, 0.2, -1.7, 0.1, -0.0, -1.2, -1.2, 0.1, -0.3, -0.5,
1954, 0.7, -0.3, 0.3, 0.6, 0.5, 0.1, 0.4, 1.1, 0.2, 0.3, 0.1, 1.4,
1955, -0.5, 1.9, 0.6, -0.1, 1.0, 1.3, 1.6, 1.5, 1.3, 1.5, 1.2, 1.0,
1956, 1.3, 1.6, 1.3, 0.9, 1.4, 1.1, 1.1, 1.2, 0.1, 1.8, 0.2, 1.1,
1957, 0.6, -0.1, 0.2, 0.2, -0.7, 0.2, 0.2, -0.5, -0.9, 0.0, -1.0, -0.3,
1958, -1.9, -0.5, 0.3, 0.4, -0.5, 0.3, 0.4, 0.9, -0.3, 0.1, -0.4, -0.6,
1959, -0.9, -1.4, 1.3, 0.4, 0.5, -0.1, -0.3, -0.1, 0.0, 0.5, 0.9, 0.9,
1960, 0.1, 0.1, 1.0, 0.8, 0.5, 0.1, 0.5, 0.8, 0.7, 0.1, 0.5, 0.8,
1961, -0.3, 0.9, -1.8, 0.8, 0.3, 0.1, 0.2, 0.2, 0.1, -0.3, 0.5, 1.5,
1962, 2.0, -0.3, 0.1, 0.2, 1.1, 0.7, 0.1, 0.6, 0.4, 1.0, 0.3, 0.2,
1963, 1.0, 0.6, 1.1, 0.8, 0.4, -0.5, -0.1, 0.0, -0.6, -1.2, -0.8, -1.2,
1964, -0.4, 0.0, 1.1, 1.1, 0.2, 0.8, 0.6, 1.5, 1.3, 1.3, 0.2, -0.3,
1965, -0.4, 0.4, 0.8, -0.5, 0.2, -0.6, -1.8, -0.7, -1.3, -0.9, -1.5, 0.2,
1966, -1.3, -0.2, -0.9, -0.2, -0.4, 0.3, 0.1, 0.6, -0.2, -0.1, -0.0, -0.3,
1967, 1.7, 1.7, 1.2, -0.0, -0.0, 0.6, 0.2, 0.7, 0.5, 0.1, -0.4, -0.6,
1968, 0.5, 1.3, 0.1, 0.0, 1.2, 1.1, 0.7, 0.3, -0.3, -0.1, -0.3, 0.2,
1969, -1.5, -0.5, 0.4, -0.4, -0.2, 0.2, -0.5, -0.1, -1.0, -0.9, -0.1, 0.4,
1970, -1.1, -1.0, 0.6, -0.1, 0.4, 1.0, -0.4, 0.6, 1.2, 1.0, 1.6, 1.9,
1971, 0.4, 2.0, 2.3, 1.7, 0.9, 0.4, 0.2, 1.5, 1.4, 1.7, 0.5, 0.3,
1972, 0.5, 1.1, 0.6, -0.1, -1.6, -0.5, -1.4, -0.5, -1.4, -0.9, -0.3, -1.3,
1973, -0.3, -1.4, 0.7, 0.1, 0.4, 1.1, 0.6, 1.3, 1.2, 0.8, 2.6, 1.8,
1974, 2.4, 2.1, 2.4, 0.9, 1.0, 0.4, 1.1, 0.8, 1.1, 0.9, -0.1, 0.2,
1975, -0.5, 0.8, 1.6, 1.2, 0.6, 1.3, 1.9, 2.0, 2.1, 1.7, 1.2, 2.1,
1976, 1.4, 1.7, 1.7, 0.3, 0.4, 0.3, -0.9, -0.8, -1.1, 0.4, 0.7, -0.3,
1977, -0.4, 1.2, -0.5, -0.4, -0.5, -0.9, -1.1, -0.8, -0.8, -1.0, -1.3, -1.1,
1978, -0.3, -2.7, -0.2, -0.3, 1.4, 0.7, 0.6, 0.4, 0.1, -0.4, -0.0, -0.1,
1979, -0.4, 1.0, 0.1, -0.1, 0.5, 0.6, 1.3, -0.2, 0.1, -0.1, -0.4, -0.7,
1980, 0.4, 0.3, -0.4, -0.6, -0.0, -0.0, -0.0, 0.4, -0.5, 0.0, -0.3, -0.1,
1981, 0.4, -0.2, -1.3, -0.1, 0.8, 1.2, 0.8, 0.7, 0.3, -0.4, 0.2, 0.5,
1982, 1.2, 0.3, 0.6, 0.1, -0.3, -1.0, -1.5, -1.7, -1.7, -1.7, -2.6, -2.2,
1983, -3.5, -3.6, -2.4, -0.9, 0.6, 0.0, -0.6, 0.1, 0.9, 0.4, -0.1, 0.0,
1984, 0.2, 0.9, -0.2, 0.3, 0.2, -0.3, 0.2, 0.4, 0.1, -0.3, 0.3, -0.1,
1985, -0.3, 1.2, 0.8, 1.2, 0.4, -0.4, -0.1, 1.0, 0.0, -0.4, -0.2, 0.2,
1986, 1.0, -1.0, 0.5, 0.3, -0.2, 1.0, 0.3, -0.4, -0.5, 0.6, -1.2, -1.4,
1987, -0.7, -1.2, -1.3, -1.4, -1.3, -1.1, -1.4, -0.9, -1.0, -0.4, -0.0, -0.5,
1988, -0.1, -0.4, 0.6, 0.1, 0.9, 0.1, 1.0, 1.5, 1.8, 1.4, 1.7, 1.2,
1989, 1.5, 1.2, 1.1, 1.6, 1.2, 0.7, 0.9, -0.3, 0.5, 0.8, -0.2, -0.5,
1990, -0.1, -1.8, -0.4, 0.2, 1.2, 0.3, 0.5, -0.2, -0.7, 0.3, -0.5, -0.2,
1991, 0.6, 0.3, -0.7, -0.6, -1.0, -0.1, 0.0, -0.4, -1.5, -1.0, -0.7, -1.8,
1992, -2.9, -0.9, -2.0, -1.0, 0.3, -0.6, -0.6, 0.4, 0.1, -1.4, -0.7, -0.6,
1993, -0.9, -0.7, -0.5, -1.2, -0.3, -0.8, -0.8, -0.9, -0.7, -1.1, -0.1, 0.2,
1994, -0.1, 0.3, -0.7, -1.3, -0.7, -0.4, -1.3, -1.2, -1.6, -1.1, -0.6, -1.2,
1995, -0.4, -0.1, 0.8, -0.7, -0.4, 0.1, 0.4, 0.3, 0.3, 0.0, 0.0, -0.5,
1996, 1.0, 0.3, 1.1, 0.8, 0.3, 1.2, 0.7, 0.7, 0.6, 0.6, -0.1, 0.9,
1997, 0.5, 1.7, -0.4, -0.6, -1.3, -1.4, -0.8, -1.4, -1.4, -1.5, -1.2, -1.0,
1998, -2.7, -2.0, -2.4, -1.4, 0.3, 1.0, 1.2, 1.2, 1.0, 1.1, 1.0, 1.4,
1999, 1.8, 1.0, 1.3, 1.4, 0.2, 0.3, 0.5, 0.4, -0.1, 1.0, 1.0, 1.4,
2000, 0.7, 1.7, 1.3, 1.2, 0.4, -0.2, -0.2, 0.7, 0.9, 1.1, 1.8, 0.8,
2001, 1.0, 1.7, 0.9, 0.2, -0.5, 0.3, -0.2, -0.4, 0.2, -0.0, 0.7, -0.8,
2002, 0.4, 1.1, -0.2, -0.1, -0.8, -0.2, -0.5, -1.0, -0.6, -0.4, -0.5, -1.1,
2003, -0.2, -0.7, -0.3, -0.1, -0.3, -0.6, 0.3, 0.1, -0.1, 0.0, -0.3, 1.1,
2004, -1.3, 1.2, 0.4, -0.9, 1.0, -0.8, -0.5, -0.3, -0.3, -0.1, -0.7, -0.8,
2005, 0.3, -3.1, 0.3, -0.6, -0.8, 0.4, 0.2, -0.3, 0.4, 1.2, -0.2, -0.0,
2006, 1.7, 0.1, 1.8, 1.1, -0.5, -0.2, -0.6, -1.0, -0.6, -1.3, 0.1, -0.3,
2007, -0.8, -0.1, 0.2, -0.1, -0.1, 0.5, -0.3, 0.4, 0.2, 0.7, 0.9, 1.7,
2008, 1.8, 2.6, 1.4, 0.7, -0.1, 0.6, 0.3, 1.0, 1.2, 1.3, 1.3, 1.4,
2009, 1.1, 1.9, 0.4, 0.8, -0.1, 0.1, 0.2, -0.2, 0.3, -1.2, -0.6, -0.7,
2010, -1.1, -1.5, -0.7, 1.2, 0.9, 0.4, 1.8, 1.8, 2.2, 1.7, 1.3, 2.9,
2011, 2.3, 2.7, 2.5, 1.9, 0.4, 0.2, 1.0, 0.4, 1.0, 0.8, 1.1, 2.5,
2012, 1.1, 0.5, 0.7, -0.3, 0.0, -0.4, -0.0, -0.2, 0.2, 0.3, 0.3, -0.6,
2013, -0.1, -0.2, 1.5, 0.2, 0.8, 1.2, 0.8, 0.2, 0.3, -0.1, 0.7, 0.1,
2014, 1.4, 0.1, -0.9, 0.8, 0.5, 0.2, -0.2, -0.7, -0.7, -0.6, -0.9, -0.6,
2015, -0.8, 0.2, -0.7, -0.0, -0.7, -0.6, -1.1, -1.4, -1.6, -1.7, -0.5, -0.6,
2016, -2.2, -2.0, -0.1, -1.2, 0.4, 0.6, 0.4, 0.7, 1.2, -0.3, -0.1, 0.3,
2017, 0.2, -0.1, 0.9, -0.2, 0.3, -0.4, 0.8, 0.5, 0.6, 0.9, 0.9, -0.1,
2018, 1.1, -0.5, 1.5, 0.5, 0.4, -0.1, 0.2, -0.3, -0.9, 0.4, -0.1, 1.0,
2019, -0.0, -1.4, -0.3, 0.1, -0.4, -0.5, -0.4, -0.1, -1.2, -0.4, -0.8, -0.6,
2020, 0.2, -0.1, -0.1, 0.2, 0.4, -0.4, 0.4, 1.1, 0.9, 0.5, 0.7, 1.8,
2021, 1.9, 1.5, 0.4, 0.3, 0.5, 0.4, 1.4, 0.6, 0.8, 0.7, 1.0, 1.5
)
#####
# Pick a date centered in each month so you can plot them as one timeline
month_centers <- data.frame(Month = toupper(month.abb),
Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
# Reformat
soi_long <- soi_std %>%
pivot_longer(names_to = "Month", values_to = "soi_z", cols = JAN:DEC) %>%
left_join(month_centers, by = "Month") %>%
mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
ribbon_fill = ifelse(soi_z > 0,"above", "below"),
ribbon_above = ifelse(soi_z > 0, soi_z, NA),
above_ymin = ifelse(soi_z > 0, 0, NA),
ribbon_below = ifelse(soi_z < 0, soi_z, NA),
below_ymin = ifelse(soi_z < 0, 0, NA)
)
# Set color scale and break points
cutpoints <- rev(seq(-4, 4, 1))
col_rdbu <- rev(RColorBrewer::brewer.pal(n = length(cutpoints), name = "RdBu"))
# fix GeomRibbon
GeomRibbon$handle_na <- function(data, params) { data }
# Make figure
ggplot(filter(soi_long, YEAR >= 1982), aes(x = Date, y = soi_z)) +
geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "La Niña"), alpha = 0.8) +
geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "El Niño"), alpha = 0.8) +
scale_fill_manual(values = c("La Niña" = col_rdbu[1], "El Niño" = col_rdbu[8])) +
geom_line(size = 0.25, color = "gray20") +
labs(x = "Date", y = "SOI", title = "ENSO Index",
caption = "Years overlapping OISST data displayed", fill = "") +
theme_gmri() +
theme(legend.position = "bottom")
When we look at average temperatures for the Gulf of Maine, and how that relates to the ENSO index, we found there was NO relationship. Meaning that a current month’s El niño/ La Nińa status did not correlate with warmer/hotter months.
North Atlantic Oscillation
A similar climate index to ENSO, with more relevance to our area is the North Atlantic Oscillation (NAO). The NAO index is based on the surface sea-level pressure difference between the Subtropical (Azores) High and the Subpolar Low.
Strong positive phases of the NAO tend to be associated with above-normal temperatures in the eastern United States and across northern Europe and below-normal temperatures in Greenland and oftentimes across southern Europe and the Middle East.
hide
# NAO Data Source:
# https://www.ncdc.noaa.gov/teleconnections/nao/
nao_idx <- tribble( #####
~"YEAR", ~"Jan", ~"Feb", ~"Mar", ~"Apr", ~"May", ~"Jun", ~"Jul", ~"Aug", ~"Sep", ~"Oct", ~"Nov", ~"Dec",
1950, 0.92, 0.40, -0.36, 0.73, -0.59, -0.06, -1.26, -0.05, 0.25, 0.85, -1.26, -1.02,
1951, 0.08, 0.70, -1.02, -0.22, -0.59, -1.64, 1.37, -0.22, -1.36, 1.87, -0.39, 1.32,
1952, 0.93, -0.83, -1.49, 1.01, -1.12, -0.40, -0.09, -0.28, -0.54, -0.73, -1.13, -0.43,
1953, 0.33, -0.49, -0.04, -1.67, -0.66, 1.09, 0.40, -0.71, -0.35, 1.32, 1.04, -0.47,
1954, 0.37, 0.74, -0.83, 1.34, -0.09, -0.25, -0.60, -1.90, -0.44, 0.60, 0.40, 0.69,
1955, -1.84, -1.12, -0.53, -0.42, -0.34, -1.10, 1.76, 1.07, 0.32, -1.47, -1.29, 0.17,
1956, -0.22, -1.12, -0.05, -1.06, 2.21, 0.10, -0.75, -1.37, 0.24, 0.88, 0.51, 0.10,
1957, 1.05, 0.11, -1.26, 0.49, -0.79, -0.72, -1.19, -0.55, -1.66, 1.32, 0.73, 0.12,
1958, -0.54, -1.06, -1.96, 0.37, -0.24, -1.38, -1.73, -1.56, -0.07, 0.16, 1.64, -0.70,
1959, -0.87, 0.68, -0.15, 0.36, 0.39, 0.40, 0.74, 0.06, 0.88, 0.89, 0.41, 0.44,
1960, -1.29, -1.89, -0.50, 1.36, 0.45, -0.21, 0.35, -1.40, 0.39, -1.73, -0.51, 0.06,
1961, 0.41, 0.45, 0.55, -1.55, -0.36, 0.86, -0.39, 0.90, 1.24, 0.51, -0.62, -1.48,
1962, 0.61, 0.55, -2.47, 0.99, -0.10, 0.16, -2.47, 0.14, -0.37, 0.41, -0.23, -1.32,
1963, -2.12, -0.96, -0.43, -1.35, 2.16, -0.43, -0.77, -0.64, 1.79, 0.94, -1.27, -1.92,
1964, -0.95, -1.43, -1.20, 0.36, 0.52, 1.29, 1.90, -1.77, 0.20, 0.74, -0.01, -0.15,
1965, -0.12, -1.55, -1.51, 0.72, -0.62, 0.29, 0.32, 0.45, 0.37, 0.38, -1.66, 1.37,
1966, -1.74, -1.39, 0.56, -0.75, 0.22, 1.05, 0.32, -1.76, -0.45, -0.68, -0.04, 0.72,
1967, -0.89, 0.19, 1.51, 0.18, -0.99, 1.40, 0.41, 1.44, 0.93, 0.07, 0.60, -0.45,
1968, 0.13, -1.29, 0.40, -1.08, -1.76, 0.33, -0.80, -0.66, -1.92, -2.30, -0.93, -1.40,
1969, -0.83, -1.55, -1.56, 1.53, 0.55, 0.55, 0.57, -1.45, 2.07, 0.66, -0.96, -0.28,
1970, -1.50, 0.64, -0.96, -1.30, 1.14, 1.55, 0.10, 0.10, -0.09, -0.92, -0.60, -1.20,
1971, -1.13, 0.24, -0.84, -0.24, 0.50, -1.57, 0.24, 1.55, 0.39, 0.58, -0.20, 0.60,
1972, 0.27, 0.32, 0.72, -0.22, 0.95, 0.88, 0.18, 1.32, -0.12, 1.09, 0.54, 0.19,
1973, 0.04, 0.85, 0.30, -0.54, -0.44, 0.39, 0.57, -0.06, -0.30, -1.24, -0.93, 0.32,
1974, 1.34, -0.14, -0.03, 0.51, -0.24, -0.14, -0.76, -0.64, 0.82, 0.49, -0.54, 1.50,
1975, 0.58, -0.62, -0.61, -1.60, -0.52, -0.84, 1.55, -0.26, 1.56, -0.54, 0.41, 0.00,
1976, -0.25, 0.93, 0.75, 0.26, 0.96, 0.80, -0.32, 1.92, -1.29, -0.08, 0.17, -1.60,
1977, -1.04, -0.49, -0.81, 0.65, -0.86, -0.57, -0.45, -0.28, 0.37, 0.52, -0.07, -1.00,
1978, 0.66, -2.20, 0.70, -1.17, 1.08, 1.38, -1.14, 0.64, 0.46, 1.93, 3.04, -1.57,
1979, -1.38, -0.67, 0.78, -1.71, -1.03, 1.60, 0.83, 0.96, 1.01, -0.30, 0.53, 1.00,
1980, -0.75, 0.05, -0.31, 1.29, -1.50, -0.37, -0.42, -2.24, 0.66, -1.77, -0.37, 0.78,
1981, 0.37, 0.92, -1.19, 0.36, 0.20, -0.45, 0.05, 0.39, -1.45, -1.35, -0.38, -0.02,
1982, -0.89, 1.15, 1.15, 0.10, -0.53, -1.63, 1.15, 0.26, 1.76, -0.74, 1.60, 1.78,
1983, 1.59, -0.53, 0.95, -0.85, -0.07, 0.99, 1.19, 1.61, -1.12, 0.65, -0.98, 0.29,
1984, 1.66, 0.72, -0.37, -0.28, 0.54, -0.42, -0.07, 1.15, 0.17, -0.07, -0.06, 0.00,
1985, -1.61, -0.49, 0.20, 0.32, -0.49, -0.80, 1.22, -0.48, -0.52, 0.90, -0.67, 0.22,
1986, 1.11, -1.00, 1.71, -0.59, 0.85, 1.22, 0.12, -1.09, -1.12, 1.55, 2.29, 0.99,
1987, -1.15, -0.73, 0.14, 2.00, 0.98, -1.82, 0.52, -0.83, -1.22, 0.14, 0.18, 0.32,
1988, 1.02, 0.76, -0.17, -1.17, 0.63, 0.88, -0.35, 0.04, -0.99, -1.08, -0.34, 0.61,
1989, 1.17, 2.00, 1.85, 0.28, 1.38, -0.27, 0.97, 0.01, 2.05, -0.03, 0.16, -1.15,
1990, 1.04, 1.41, 1.46, 2.00, -1.53, -0.02, 0.53, 0.97, 1.06, 0.23, -0.24, 0.22,
1991, 0.86, 1.04, -0.20, 0.29, 0.08, -0.82, -0.49, 1.23, 0.48, -0.19, 0.48, 0.46,
1992, -0.13, 1.07, 0.87, 1.86, 2.63, 0.20, 0.16, 0.85, -0.44, -1.76, 1.19, 0.47,
1993, 1.60, 0.50, 0.67, 0.97, -0.78, -0.59, -3.18, 0.12, -0.57, -0.71, 2.56, 1.56,
1994, 1.04, 0.46, 1.26, 1.14, -0.57, 1.52, 1.31, 0.38, -1.32, -0.97, 0.64, 2.02,
1995, 0.93, 1.14, 1.25, -0.85, -1.49, 0.13, -0.22, 0.69, 0.31, 0.19, -1.38, -1.67,
1996, -0.12, -0.07, -0.24, -0.17, -1.06, 0.56, 0.67, 1.02, -0.86, -0.33, -0.56, -1.41,
1997, -0.49, 1.70, 1.46, -1.02, -0.28, -1.47, 0.34, 0.83, 0.61, -1.70, -0.90, -0.96,
1998, 0.39, -0.11, 0.87, -0.68, -1.32, -2.72, -0.48, -0.02, -2.00, -0.29, -0.28, 0.87,
1999, 0.77, 0.29, 0.23, -0.95, 0.92, 1.12, -0.90, 0.39, 0.36, 0.20, 0.65, 1.61,
2000, 0.60, 1.70, 0.77, -0.03, 1.58, -0.03, -1.03, -0.29, -0.21, 0.92, -0.92, -0.58,
2001, 0.25, 0.45, -1.26, 0.00, -0.02, -0.20, -0.25, -0.07, -0.65, -0.24, 0.63, -0.83,
2002, 0.44, 1.10, 0.69, 1.18, -0.22, 0.38, 0.62, 0.38, -0.70, -2.28, -0.18, -0.94,
2003, 0.16, 0.62, 0.32, -0.18, 0.01, -0.07, 0.13, -0.07, 0.01, -1.26, 0.86, 0.64,
2004, -0.29, -0.14, 1.02, 1.15, 0.19, -0.89, 1.13, -0.48, 0.38, -1.10, 0.73, 1.21,
2005, 1.52, -0.06, -1.83, -0.30, -1.25, -0.05, -0.51, 0.37, 0.63, -0.98, -0.31, -0.44,
2006, 1.27, -0.51, -1.28, 1.24, -1.14, 0.84, 0.90, -1.73, -1.62, -2.24, 0.44, 1.34,
2007, 0.22, -0.47, 1.44, 0.17, 0.66, -1.31, -0.58, -0.14, 0.72, 0.45, 0.58, 0.34,
2008, 0.89, 0.73, 0.08, -1.07, -1.73, -1.39, -1.27, -1.16, 1.02, -0.04, -0.32, -0.28,
2009, -0.01, 0.06, 0.57, -0.20, 1.68, -1.21, -2.15, -0.19, 1.51, -1.03, -0.02, -1.93,
2010, -1.11, -1.98, -0.88, -0.72, -1.49, -0.82, -0.42, -1.22, -0.79, -0.93, -1.62, -1.85,
2011, -0.88, 0.70, 0.61, 2.48, -0.06, -1.28, -1.51, -1.35, 0.54, 0.39, 1.36, 2.52,
2012, 1.17, 0.42, 1.27, 0.47, -0.91, -2.53, -1.32, -0.98, -0.59, -2.06, -0.58, 0.17,
2013, 0.35, -0.45, -1.61, 0.69, 0.57, 0.52, 0.67, 0.97, 0.24, -1.28, 0.90, 0.95,
2014, 0.29, 1.34, 0.80, 0.31, -0.92, -0.97, 0.18, -1.68, 1.62, -1.27, 0.68, 1.86,
2015, 1.79, 1.32, 1.45, 0.73, 0.15, -0.07, -3.18, -0.76, -0.65, 0.44, 1.74, 2.24,
2016, 0.12, 1.58, 0.73, 0.38, -0.77, -0.43, -1.76, -1.65, 0.61, 0.41, -0.16, 0.48,
2017, 0.48, 1.00, 0.74, 1.73, -1.91, 0.05, 1.26, -1.10, -0.61, 0.19, -0.00, 0.88,
2018, 1.44, 1.58, -0.93, 1.24, 2.12, 1.09, 1.39, 1.97, 1.67, 0.93, -0.11, 0.61,
2019, 0.59, 0.29, 1.23, 0.47, -2.62, -1.09, -1.43, -1.17, -0.16, -1.41, 0.28, 1.20,
2020, 1.34, 1.26, 1.01, -1.02, -0.41, -0.15, -1.23, 0.12, 0.98, -0.65, 2.54, -0.30,
2021, -1.11, 0.14, 0.73, -1.43, -1.24, 0.77, 0.03, -0.28, -0.21, -2.29, -0.18, 0.29
)
#####
# Reformat and add a center date
month_centers <- data.frame(Month = month.abb,
Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
# Reformat
nao_long <- nao_idx %>%
pivot_longer(names_to = "Month", values_to = "nao_z", cols = Jan:Dec) %>%
left_join(month_centers, by = "Month") %>%
mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
ribbon_fill = ifelse(nao_z > 0,"above", "below"),
ribbon_above = ifelse(nao_z > 0, nao_z, NA),
above_ymin = ifelse(nao_z > 0, 0, NA),
ribbon_below = ifelse(nao_z < 0, nao_z, NA),
below_ymin = ifelse(nao_z < 0, 0, NA)
)
# Make figure
ggplot(filter(nao_long, YEAR >= 1982),
aes(x = Date, y = nao_z)) +
geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "Positive Forcing"), alpha = 0.8) +
geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "Negative Forcing"), alpha = 0.8) +
scale_fill_manual(values = c("Positive Forcing" = col_rdbu[1], "Negative Forcing" = col_rdbu[8])) +
geom_line(size = 0.25, color = "gray20") +
labs(x = "Date", y = "NAO", title = "NAO Index",
caption = "Years overlapping OISST data displayed", fill = "") +
theme_gmri() +
theme(legend.position = "bottom")
Looking at the NAO index, we found there was also not a strong relationship to our region’s Sea surface temperature anomalies. This lends strength to the argument that recent changes are independent of normal climate oscillations, but is not the end of the story.
Daily Temperature Record
The Gulf of Maine has a regular seasonal cycle hitting its lowest temperatures in March, and experiencing its peak highs in August. The difference between these two seasons is ~\(20-25^oF\) or ~\(11-12^oC\) (eyeballing it here).
hide
# Last Calendar Year
this_yr <- region_hw %>%
filter(year(time) == 2021)
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"MHW Threshold" = "coral3",
"Daily Climatology" = "gray30")
# Set the label with degree symbol
ylab <- "Sea Surface Temperature"
# Plot the last 365 days
hw_temp_p <- ggplot(this_yr, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
# geom_textpath(aes(y = mhw_thresh, color = "MHW Threshold"), label = "Heatwave Threshold", hjust = 0.8, lty = 3) +
#geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
geom_textpath(aes(y = seas, color = "Daily Climatology"), label = "Climatology Mean", hjust = 0.5, lty = 2) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
theme(legend.title = element_blank(),
legend.position = "none") +
labs(x = "",
y = ylab)
# Plot the last 365 days - anomaly scale
linetype_key <- c(
"Sea Surface Temperature Anomaly" = 1,
"Heatwave Event" = 1,
"MHW Threshold" = 3,
"Daily Climatology" = 2)
# Same Plot as Anomalies
hw_anom_p <- this_yr %>%
mutate(
sst = sst,
seas = seas,
sst_anom = sst_anom,
mhw_thresh = mhw_thresh,
anom_thresh = mhw_thresh - seas,
anom_hwe = hwe - seas) %>%
ggplot(aes(x = time)) +
geom_segment(aes(x = time, xend = time,
y = 0, yend = sst_anom),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time,
y = anom_thresh, yend = anom_hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
# geom_textpath(aes(y = 0, color = "Daily Climatology"),
# label = "Climatology Mean", hjust = 0.5, lty = 2, show.legend = FALSE) +
scale_color_manual(values = color_vals) +
scale_linetype_manual(values = linetype_key, guide = "none") +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
theme(legend.title = element_blank(),
legend.position = "top") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
labs(x = "",
y = "Temperature Anomaly",
caption = paste0("Climate reference period : 1982-2011"))
# Show figure
hw_temp_p / hw_anom_p
2021 Heatwave Events
hide
# Get Event numbers
events_21 <- this_yr %>% drop_na(mhw_event_no) %>% distinct(mhw_event_no) %>% pull(mhw_event_no)
events_num <- length(events_21)
# Pull their data
events_dat <- filter(region_hw, mhw_event_no %in% as.character(events_21),
yr <= 2021)
# Days in HW state
days_hw <- nrow(events_dat)
# How long on average
avg_duration <- round(days_hw / events_num, digits = 0)
During 2021 there were 5 main marine heatwave events, including one that carried over from the previous year and one that has continued into 2022. The average duration for these events was 72 days (excludes ongoing days in 2022).
Heatwave Seasonality
Over the last ~40 years the frequency, duration, and intensity of heatwave events has increased. This year was no exception with 358 days meeting the criteria for heatwave status.
hide
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# Polar plot comparing when each year experienced peak heatwave conditions
polar_dat <- grid_data %>% filter(year %in% c("2012", "2021")) %>%
mutate(hw_line = ifelse(mhw_event == TRUE, sst_anom, NA),
hw_mag = ifelse(year == "2012", hw_line * -1, hw_line))
hide
# Things to highlight:
# when the heatwaves ocurred, how they overlapped
# severity?
# Side by side - need data separately
polar_dat %>%
ggplot(aes(x = flat_date, xend = flat_date,
y = sst_anom, yend = 0, color = mhw_event)) +
geom_segment(alpha = 0.9, show.legend = F) +
geom_textpath(aes(x = flat_date, y = mhw_thresh - seas),
label = "Heatwave Threshold", linetype = 1, color = "black", hjust = 0.85, straight = T) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(mult = c(0,0))) +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C"),
expand = expansion(mult = c(0,0))) +
facet_wrap( ~ year, nrow = 2) +
scale_color_manual(values = c("TRUE" = "darkred", "FALSE" = "black")) +
# scale_color_gmri() +
# theme_gmri() +
labs(x = "", y = "Temperature Anomalies", fill = "Year")
Following Each Heatwave
Do some breakdowns of how each heatwave looked: - where was heat concentrated - how long did it last - was anything unique about it?
Heatwave Frequency Heatmap
hide
# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) * c(-1, 1)
# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
labs(x = "",
y = "",
"\nClimate reference period : 1982-2011") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent",
limit = limit) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
# Assemble pieces
heatwave_heatmap
Temperature Horizons
hide
# make slight modification to grid_data
hori_data <- grid_data %>%
mutate(year = factor(year),
year = fct_rev(year))
# cut out outliers and set cutpoints
cutpoints <- hori_data %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# Set origin
ori <- 0
# Manually pick cut points
sca <- seq(-4, 4, length.out = 9)
sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "\u00b0C")
sca_labels <- rev(sca_bins)
# Function to plot one or more years
anom_horizon_plot <- function(hori_data, origin = ori, scale_cutpoints = sca, labels = sca_labels){
ggplot(hori_data) +
geom_horizon(aes(flat_date,
sst_anom,
fill = ..Cutpoints..),
origin = ori,
horizonscale = sca) +
scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
facet_grid(year~.) +
theme(
panel.grid = element_blank(),
panel.spacing.y=unit(0, "lines"),
strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
legend.position = "left"
) +
scale_x_date(expand=c(0,0),
date_breaks = "1 month",
date_labels = "%b") +
labs(x = '',
fill = "Temperature Anomaly \u00b0C")
}
# Plot Horizon Plot Using All Years
anom_horizon_plot(hori_data = filter(hori_data, yr <= 2021),
origin = ori,
scale_cutpoints = sca) +
labs(title = "Temperature Anomaly Horizons",
subtitle = "How deviation from the norm has become the new-normal",
caption = "Climatatological Reference Period: 1982-2011")
Comparing the Gulf’s Two Hottest Years
Early into 2021 it was apparent that the year was on-par with the previous title-holder for warmest year on record.
hide
# Reset the outliers
# cut out outliers and set cutpoints
cutpoints_2 <- hori_data %>%
filter(year %in% c("2012", "2021")) %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# 2012
hori_2012 <- hori_data %>%
filter(year == "2012") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top")
# 2021
hori_2021 <- hori_data %>% filter(year == "2021") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top") +
theme(plot.caption = element_text(colour = "gray25", size = 6)) +
labs(caption = "Climatatological Reference Period: 1982-2011")
# Put them together as one figure
comparison_horizons <- hori_2012 / hori_2021 + plot_layout(guides = "collect")
# Do a bunch of formatting
comparison_horizons &
theme(
#Legend
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines")) &
guides(fill = guide_legend(title = "Temperature Anomaly",
title.hjust = 0.5,
title.position = "top",
label.position = "bottom",
nrow = 1,
byrow = T, reverse = T)) &
plot_annotation(title = "Comparing the Gulf of Maine's Two Hottest Years",
subtitle = "Duration of temperature anomaly horizons, colored by their strength")
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Streamplot
hide
# Set Factor Order for Bins
bin_order <- c(
"Less than -1",
"-1 to -0.5",
"Within 0.5",
"0.5 to 1",
"Greater than 1")
# Set Colors
bin_colors <- RColorBrewer::brewer.pal(n = length(bin_order),
name = "RdBu")
bin_colors <- rev(bin_colors)
# Count Days in Bins
day_types <- region_hw %>%
mutate(
`Less than -1` = if_else(sst_anom <= -1, TRUE, FALSE),
`-1 to -0.5` = if_else(between(sst_anom, -1, -0.5), TRUE, FALSE),
`Within 0.5` = if_else(between(sst_anom, -0.5, 0.5), TRUE, FALSE),
`0.5 to 1` = if_else(between(sst_anom, 0.5, 1), TRUE, FALSE),
`Greater than 1` = if_else(sst_anom >= 1, TRUE, FALSE)) %>%
group_by(yr) %>%
summarise(
`Less than -1` = sum(`Less than -1`, na.rm = T),
`-1 to -0.5` = sum(`-1 to -0.5`, na.rm = T),
`Within 0.5` = sum(`Within 0.5`, na.rm = T),
`0.5 to 1` = sum(`0.5 to 1`, na.rm = T),
`Greater than 1` = sum(`Greater than 1`, na.rm = T),
.groups = "drop") %>%
pivot_longer(cols = "Less than -1":"Greater than 1", names_to = "Anomaly Strength", values_to = "day_total", values_drop_na = FALSE) %>%
mutate(day_total = if_else(is.na(day_total), 0L, day_total),
`Anomaly Strength` = factor(`Anomaly Strength`, levels = bin_order))
# Plot the Streamplot
hw_streamplot <- ggplot(day_types, aes(yr, y = day_total, fill = `Anomaly Strength`)) +
geom_stream(type = "proportion") +
#geom_segment(y = .5, yend = 0.5, linetype = 3, color = "gray50", x = 1981, xend = 2021, size = 0.25) +
scale_x_continuous(breaks = seq(1980,2030, by = 10),
minor_breaks = seq(1975, 2030, by = 5),
expand = c(.03,0)) +
scale_y_continuous(labels = scales::percent_format(), expand = c(0,0)) +
scale_fill_manual(values = setNames(bin_colors, bin_order)) +
labs(y = "", x = "", fill = "",
title = "A Warmer Gulf of Maine",
subtitle = "Fractions of Each Year Experienced as Hot/Cold Anomalies",
caption = "1981 not a complete years in the data*") +
theme_minimal(base_size = 10) +
theme(
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines"),
panel.grid.minor = element_blank(),
panel.grid = element_line(size = 0.2),
plot.margin = unit(c(.5,1,.3,.5), "cm"),
plot.title.position = "plot",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
legend.margin = margin(c(7,0,0,0)),
legend.justification = "center") +
guides(fill = guide_legend(title = "Anomaly Strength \u00b0C",
title.hjust = 0.5,
title.position = "left",
label.position = "bottom",
nrow = 1,
byrow = T))
# Show figure
hw_streamplot
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Characteristics of the Hottest Years
Multivariate figure comparing characteristics of both years: 2012 & 2021
hide
# Use all years so they can be ranked independently
# Metrics?
year_mets <- hori_data %>%
filter(year != "2022") %>%
group_by(year) %>%
summarise(
total_days = n(),
`Average Temperature` = round(mean(sst, na.rm = T), 2),
`Temperature Above Average` = round(mean(sst_anom, na.rm = T), 2),
`Cumulative Degrees Above Average` = round(sum(sst_anom, na.rm = T), 0),
hw_days = sum(mhw_event, na.rm = T),
hw_events = n_distinct(mhw_event_no),
peak_temp = max(sst),
peak_anom = max(sst_anom)) %>%
mutate(
hw_day_pct = round((hw_days/total_days) * 100, 2),
`Average Heatwave Length` = round(hw_days/hw_events, 1)
) %>%
pivot_longer(names_to = "Var", values_to = "Metric", cols = `Average Temperature`:`Average Heatwave Length`)
# Snipe out the top5 for each
year_met_ranks <- year_mets %>%
filter(Var %in% c("Average Temperature", "Temperature Above Average", "Average Heatwave Length", "Cumulative Degrees Above Average")) %>%
group_by(Var) %>%
slice_max(n = 5, order_by = Metric) %>%
ungroup() %>%
mutate(
yr_focus = ifelse(year == "2021", TRUE, FALSE),
year = reorder_within(year, Metric, Var))
year_met_ranks %>%
ggplot(aes(year, Metric, color = yr_focus)) +
geom_col(aes(fill = yr_focus), show.legend = F) +
geom_label(aes(label = Metric), show.legend = FALSE) +
facet_wrap(~Var, scales = "free", ncol = 2) +
coord_flip() +
scale_x_reordered() +
scale_y_continuous(expand = expansion(mult = c(0, .2))) +
labs(y = "",
x = NULL,
title = "How 2021 Stacks Up:",
subtitle = "Ranking Characteristics of Acute and Prolonged Warming Events") +
theme(plot.title = element_text(hjust = 0)) +
scale_color_gmri() +
scale_fill_gmri() +
theme(panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major. = element_line(linetype = 3, color = "gray50"))
Mapping Record Temperatures
2021 was an exceptionally warm year for many parts of the World, not just the Gulf of Maine. By taking the average temperature for the year for all years, we can rank how this year compares to previous years in the OISST data. For many places around the world this was one of the top 3 warmest years on record.
hide
# Load the rankings done in python, do some checking
ranks_path <- str_c(oisst_path, 'temp_rankings/yrly_ranks_1982to2021.nc')
# Pull this year
ranks_21 <- stack(ranks_path)[["X2021"]]
# Reclassify Ranks
# Ranks are all integers so we just need to be sure that they are captured in the bounds we want
reclass_df <- c(
0, 1.1, -1, # Coldest Year
1.1, 2.1, -2, # Second Coldest
2.1, 3.1, -3, # Third Coldest
3.1, 4.1, -4, # Fourth Coldest
4.1, 5.1, -5, # Fifth Coldest
3.1, 35, NA, # Things to Mask
35.1, 36, 5, # Fifth Warmest
36.1, 37, 4, # Fourth Warmest
37.1, 38, 3, # Third Warmest
38.1, 39, 2, # Second Warmest
39.1, 40, 1 # Warmest
)
# make reclassification matrix
reclass_m <- matrix(reclass_df, ncol = 3, byrow = TRUE)
# reclassify - masking middle values handled with reclassification
ranks_reclass <- reclassify(ranks_21, reclass_m)
# Make Stars for ggplot
ranks_st <- st_as_stars(rotate(ranks_reclass))
# Mutate to get levels as factors for legend
ranks_st <- ranks_st %>%
mutate(
temp_rank = case_when(
X2021 == -1 ~ "Coldest",
X2021 == -2 ~ "2nd Coldest",
X2021 == -3 ~ "3rd Coldest",
X2021 == -4 ~ "4th Coldest",
X2021 == -5 ~ "5th Coldest",
X2021 == 5 ~ "5th Warmest",
X2021 == 4 ~ "4th Warmest",
X2021 == 3 ~ "3rd Warmest",
X2021 == 2 ~ "2nd Warmest",
X2021 == 1 ~ "Warmest",
TRUE ~ "Not a Record Year"),
temp_rank = fct_drop(temp_rank),
temp_rank = factor(temp_rank, levels = c(
"Warmest", "2nd Warmest", "3rd Warmest", "4th Warmest", "5th Warmest",
"Not a Record Year",
"5th Coldest", "4th Coldest", "3rd Coldest", "2nd Coldest", "Coldest"))) %>%
select(-X2021)
# Make Map
ranks_map <- ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world_sf, size = .25) +
scale_fill_brewer(palette = "RdBu", na.value = "transparent") +
scale_x_continuous(breaks = seq(-180, 180, 90),
expand = expansion(mult = c(0,0))) +
scale_y_continuous(breaks = seq(-90, 90, 30),
expand = expansion(mult = c(0,0))) +
labs(title = "Record Setting Temperatures of 2021",
subtitle = "Regions that experienced record setting temperatures",
fill = "",
caption = "Temperature rankings cover period of 1982-2021") +
map_theme
ranks_map
hide
# Zoom in on Northeast US
new_england <- ne_states("united states of america", returnclass = "sf")
canada <- ne_states("canada", returnclass = "sf")
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = new_england, size = .25) +
geom_sf(data = canada, size = .25) +
scale_fill_brewer(palette = "RdBu", na.value = "transparent")+
coord_sf(xlim = c(-80, -55), ylim = c(34, 49)) +
labs(x = "", y = "", fill = "2021 Temperature Rank") +
map_theme
hide
# Mapping in orthographic projection:
# Function for making a bounding box
make_bbox <- function(xlim = c(-75, -60), ylim = c(32, 46), crs = 4326){
# Build bbox to crop with
crop_bbox <- st_bbox(
c(xmin = xlim[1],
ymin = ylim[1],
xmax = xlim[2],
ymax = ylim[2]),
crs = st_crs(crs))
return(crop_bbox)
}
#zoom_bbox in different CRS
make_bbox() %>%
st_transform(ortho)
# # Issues: The warping step is having issues, try cropping first
# crop_grid_bydim <- function(in_ras, xlim = c(-80, -50), ylim = c(20, 60), crs = crs){
#
#
#
# st_crop(in_ras, crop_bbox)
# }
#
#
# # Warp it to orthographic projection
# warp_to_ortho <- function(in_grid_st, ortho_lat = 25, ortho_lon = -50){
# # orthographic projection
# lat <- ortho_lon
# lon <- ortho_lon
# ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')
#
# # Build grid in the crs you wish to warp to
# projection_grid <- in_grid_st %>%
# st_transform(crs = ortho) %>%
# st_bbox() %>%
# st_as_stars()
#
# # Warp to projection grid of chosen CRS
# region_warp_ras <- in_grid_st %>%
# st_warp(projection_grid)
#
# return(region_warp_ras)
#
# }
#
#
# # Crop it down to a manageable area
# ranks_cropped <- crop_grid_bydim(in_ras = st_as_stars(rotate(ranks_21)), crs = 4326)
#
# # use warping function:
# ranks_ortho <- warp_to_ortho(in_grid_st = ranks_cropped)
# source:
# https://gist.github.com/rafapereirabr/26965dd851debad32ad2e659024ba451
# Orthographic projection
ortho_lat <- 25
ortho_lon <- -60
ortho <- paste0('+proj=ortho +lat_0=', ortho_lat, ' +lon_0=', ortho_lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')
# Can we just project the original?
ranks_ortho <- projectRaster(rotate(ranks_21), crs = ortho)
ortho_st <- st_as_stars(ranks_ortho) %>%
mutate(
temp_rank = case_when(
X2021 == 1 ~ "Coldest",
X2021 == 2 ~ "2nd Coldest",
X2021 == 3 ~ "3rd Coldest",
X2021 == 4 ~ "4th Coldest",
X2021 == 5 ~ "5th Coldest",
X2021 == 36 ~ "5th Warmest",
X2021 == 37 ~ "4th Warmest",
X2021 == 38 ~ "3rd Warmest",
X2021 == 39 ~ "2nd Warmest",
X2021 == 40 ~ "Warmest",
TRUE ~ "Not a Record Year"),
temp_rank = fct_drop(temp_rank),
temp_rank = factor(temp_rank, levels = c(
"Warmest", "2nd Warmest", "3rd Warmest", "4th Warmest", "5th Warmest",
"Not a Record Year",
"5th Coldest", "4th Coldest", "3rd Coldest", "2nd Coldest", "Coldest")))
# Try to:
# Fix polygons so they don't get cut in ortho projection
world <- rnaturalearth::ne_countries(scale = 'small', returnclass = 'sf')
world_ortho <- st_cast(world, 'MULTILINESTRING') %>%
st_cast('LINESTRING', do_split=TRUE) %>%
mutate(npts = mapview::npts(geometry, by_feature = TRUE)) %>%
st_cast('POLYGON') %>%
st_transform(ortho) %>%
filter( is.na(st_dimension(.)) == FALSE )
# and map
ggplot() +
geom_stars(data = select(ortho_st, -X2021)) +
scale_fill_brewer(palette = "RdBu", na.value = "gray50") +
coord_sf(crs = ortho, expand = F) +
map_theme
hide
# Preserving polygons in ortho view:
# https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1
# checking geometries https://r-spatial.org/r/2017/03/19/invalid.html#empty-geometries
# Define the orthographic projection
# Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180
lat <- 45
lon <- -65
ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs')
# Load the world outlines
world <- ne_states(returnclass = "sf")
# Define the polygon that will help you finding the "blade"
# to split what lies within and without your projection
circle <- st_point(x = c(0,0)) %>%
st_buffer(dist = 6371000) %>%
st_sfc(crs = ortho)
# Project this polygon in lat-lon
circle_longlat <- circle %>% st_transform(crs = 4326)
# circle_longlat cannot be used as it is
# You must decompose it into a string with ordered longitudes
# Then complete the polygon definition to cover the hemisphere
if(lat != 0) {
circle_longlat <- st_boundary(circle_longlat)
circle_coords <- st_coordinates(circle_longlat)[, c(1,2)]
circle_coords <- circle_coords[order(circle_coords[, 1]),]
circle_coords <- circle_coords[!duplicated(circle_coords),]
# Rebuild line
circle_longlat <- st_linestring(circle_coords) %>% st_sfc(crs = 4326)
if(lat > 0) {
rectangle <- list(rbind(circle_coords,
c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
c(X = 180, Y = 90),
c(X = -180, Y = 90),
c(X = -180, circle_coords[1, 'Y']),
circle_coords[1, c('X','Y')])) %>%
st_polygon() %>% st_sfc(crs = 4326)
} else {
rectangle <- list(rbind(circle_coords,
c(X = 180, circle_coords[nrow(circle_coords), 'Y']),
c(X = 180, Y = -90),
c(X = -180, Y = -90),
c(X = -180, circle_coords[1, 'Y']),
circle_coords[1, c('X','Y')])) %>%
st_polygon() %>% st_sfc(crs = 4326)
}
circle_longlat <- st_union(st_make_valid(circle_longlat), st_make_valid(rectangle))
}
# This visualization shows the visible emisphere in red
ggplot() +
geom_sf(data = world) +
geom_sf(data = circle_longlat, color = 'red', fill = 'red', alpha = 0.3)
# A small negative buffer is necessary to avoid polygons still disappearing in a few pathological cases
# I should not change the shapes too much
visible <- st_intersection(st_make_valid(world),
st_buffer(circle_longlat, -0.09)) %>%
st_transform(crs = ortho)
# DISCLAIMER: This section is the outcome of trial-and-error and I don't claim it is the best approach
# Resulting polygons are often broken and they need to be fixed
# Get reason why they're broken
broken_reason <- st_is_valid(visible, reason = TRUE)
# First fix NA's by decomposing them
# Remove them from visible for now
na_visible <- visible[is.na(broken_reason),]
visible <- visible[!is.na(broken_reason),]
# Open and close polygons
na_visible <- st_cast(na_visible, 'MULTILINESTRING') %>%
st_cast('LINESTRING', do_split=TRUE)
na_visible <- na_visible %>%
mutate(npts = mapview::npts(geometry, by_feature = TRUE))
# Exclude polygons with less than 4 points
na_visible <- na_visible %>%
filter(npts >=4) %>%
select(-npts) %>%
st_cast('POLYGON')
# Fix other broken polygons
broken <- which(!st_is_valid(visible))
for(land in broken) {
result = tryCatch({
# visible[land,] <- st_buffer(visible[land,], 0) # Sometimes useful sometimes not
visible[land,] <- st_make_valid(visible[land,]) %>%
st_collection_extract()
}, error = function(e) {
visible[land,] <<- st_buffer(visible[land,], 0)
})
}
# Bind together the two tables
visible <- rbind(visible, na_visible)
# Final plot
ggplot() +
geom_sf(data = circle,
#fill = 'aliceblue') + # if you like the color
fill = NA) +
geom_sf(data = st_collection_extract(visible)) +
coord_sf(crs = ortho)
Comparing Regions
We are often asked when speaking to the degree of the Gulf of Maine is warming: > “Where is it warming faster?”
This question is challenging because the Gulf of Maine, and any area in the ocean has a unique size/shape/dynamic that makes them difficult to directly compare directly. To put the Gulf of Maine’s warming into perspective against similar ecologically and oceanographically relevant units we can compare it against the large marine ecosystems (LME) of the world.
hide
# Load the LME timeseries to calculate rates of warming
lme_names <- get_region_names("lme") # names
lme_paths <- get_timeseries_paths("lme", mac_os = "mojave") # paths
# Data
lme_oisst <- map(lme_names, ~ read_csv(lme_paths[[.x]][["timeseries_path"]], col_types = cols(), guess_max = 1e6))
lme_oisst <- setNames(lme_oisst, lme_names)
# Add Gulf of Maine Timeseries to the List with LME's
lme_oisst[["gulf_of_maine"]] <- region_timeseries
# Make sure dates are formatted
lme_oisst <- map(lme_oisst, ~ mutate(.x, time = as.Date(time)))
#### Warming Rates ####
lme_8221 <- map(lme_oisst, ~ mutate(.x, year = year(time)) %>% filter(year %in% c(1982:2021)))
# Get warming rates
lme_rates <- map_dfr(lme_8221, function(lme_data){
# Uses area weighted SST
lme_yearly <- group_by(lme_data, year) %>%
summarise(sst = mean(area_wtd_sst, na.rm = T))
yrly_warming <- lm(sst ~ year, data = lme_yearly)
mod_details <- broom::tidy(yrly_warming)
yrly_rate <- mod_details[2,"estimate"]
names(yrly_rate) <- "annual_warming_rate_C"
yrly_rate
}, .id = "Region") %>%
arrange(desc(annual_warming_rate_C))
# Pull out the top 13
rates_ordered <- lme_rates %>%
mutate(`Warming Rank` = row_number(),
`Warming Rate F` = as_fahrenheit(annual_warming_rate_C, data_type = "anomalies"),
Region = str_replace_all(Region, "_", " "),
Region = str_to_title(Region)) %>%
select(Region, `Warming Rank`, `Warming Rate C` = annual_warming_rate_C, `Warming Rate F`)
# symbol for degree: \u00b0
hide
# Make table of top 10*
rates_ordered %>%
slice(1:10) %>%
mutate_if(is_numeric, round, 3) %>%
gt(rowname_col = "Region") %>%
tab_header(
title = md("**Warming Rates of Large Marine Ecosystems**")) %>%
tab_stubhead(label = "Region") %>%
tab_style(
style = list(
cell_fill(color = "gray90"),
cell_text(style = "oblique")),
locations = cells_body(
rows = Region == "Gulf Of Maine")) %>%
tab_spanner(label = md("*Annual Temperature Change*"),
columns = c(`Warming Rate C`, `Warming Rate F`)) %>%
cols_label(
`Warming Rate F` = "\u00b0F",
`Warming Rate C` = "\u00b0C") %>%
tab_source_note(source_note = md("**Data Source:** *NOAA OISSTv2 Daily Sea Surface Temperature Data.*")) %>%
tab_source_note(source_note = md("**Notes:** *Temperatures adjusted for latitudinal distotion.*") ) %>%
tab_footnote(
footnote = "Gulf of Maine not a true large marine ecosystem",
locations = cells_stub(rows = which(rates_ordered$Region == "Gulf Of Maine"))
)
| Warming Rates of Large Marine Ecosystems | |||
|---|---|---|---|
| Region | Warming Rank | Annual Temperature Change | |
| °C | °F | ||
| Baltic Sea | 1 | 0.047 | 0.085 |
| Gulf Of Maine1 | 2 | 0.044 | 0.080 |
| Black Sea | 3 | 0.044 | 0.079 |
| Scotian Shelf | 4 | 0.041 | 0.073 |
| Iceland Shelf And Sea | 5 | 0.038 | 0.069 |
| Northeast Us Continental Shelf | 6 | 0.038 | 0.068 |
| North Sea | 7 | 0.035 | 0.063 |
| Norwegian Sea | 8 | 0.032 | 0.057 |
| Sea Of Japan | 9 | 0.031 | 0.056 |
| Mediterranean Sea | 10 | 0.030 | 0.055 |
| Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data. | |||
| Notes: Temperatures adjusted for latitudinal distotion. | |||
|
1
Gulf of Maine not a true large marine ecosystem
|
|||
hide
#
rates_path <- str_c("{box_root}RES_Data/OISST/oisst_mainstays/warming_rates/annual_warming_rates{trend_start_year}to{trend_end_year}.nc")
# Get paths to the LME's and their timeseries
lme_paths <- gmRi::get_timeseries_paths(region_group = "lme", mac_os = "mojave")
Biological Events
- basking sharks
- hammerheads
- blue crab
A work by Adam A. Kemberling
Akemberling@gmri.org